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ABSTRACT 


Research  results  summarizing  the  past  two  years’  efforts  are  presented.  These  include 
analysis  of  both  the  large-  and  small-scale  equations  of  additive  turbulent  decomposition  (ATD) 
for  the  2-D  incompressible  Navier-Stokes  equations.  Both  of  these  solution  procedures  are 
complete  for  Cartesian  coordinates  and  check  out  is  essentially  fmished  for  generalized 
coordinates.  We  discuss  our  method  used  to  filter  solutions  to  the  large-scale  equations,  and  the 
domain  decomposition  methods  we  are  employing  to  recouple  local  small-scale  solutions  to 
produce  global  solutions  on  the  small  scale.  We  also  present  results  obtained  from  our  large- 
scale  equations/chaotic  map  turbulence  modeling  approach  for  intermittent  flow. 


Aooesslon  For 

RTIS  GRA&l 

ca^ 

DTIC  TAB 

□ 

Unaiuiouncod 

Jma  t  i  j-'  1  vi  at  i  on _ 

□ 

By - 

Dlatrlbiitiop/' 


ATailnfclilty  Oodea 
|AT»ii  aad/or 
Special 


2 


1.  INTRODUCTION 


This  report  surnmarizes  research  carried  out  under  AFOSR  Grant  #F49620-92-J-0113 
“Additive  Turbulent  Decomposition  of  the  Incompressible  and  Compressible  Navier-Stokes 
Equations,”  during  the  period  15  Nov  1991  to  14  Nov  1993.  The  work  originally  proposed  for 
this  grant  included  continuing  studies  of  additive  turbulent  decomposition  (ATD)  in  the  context 
of  two-dimensional  (2-D)  incompressible  flows  in  generalized  coordinates,  and  initial  studies  of 
2-D  compressible  turbulent  flows.  Due  to  a  significant  reduction  in  funding  from  the  proposed 
amount,  the  second  of  these  tasks  has  not  been  carried  out  under  this  grant  We  note,  however, 
that  this  work  has  begun  under  a  separate  AFOSR  Grant  #F49620-92-J-0488,  effective  1  Sep 
1992,  and  will  be  reported  at  the  appropriate  time. 

Studies  completed  during  the  first  year  of  this  project  have  been  discussed  in  fair  detail  in 
the  Annual  Report  for  1991-1992  (McDonough  et  al.,  1992),  so  this  work  will  mainly  only  be 
summarized  in  the  present  report  for  the  sake  of  completeness.  This  research  included  the 
analyses  required  for  the  construction  of  ATD  in  generalized  coordinates  and  discussion  of  the 
projection  method  algorithm  employed  to  solve  the  large-scale  equations.  Details  of  the  small- 
scale  equations  and  their  solutions  were  presented,  and  the  complete  ATD  solution  algorithm 
was  discussed  briefly.  Two  key  parts  of  this  were  treatment  of  aliasing  of  the  coarsely  resolved 
solutions  of  the  large-scale  equations;  and  parallelization  of  ATD.  We  will  include  considerably 
more  regarding  both  of  these  topics  in  the  present  report,  as  well  as  on  new  domain 
decomposition  approaches  for  coupling  the  local  small-scale  solutions. 

A  second  major  topic  was  also  presented  in  the  referenced  Annual  Report,  namely 
turbulence  models  based  on  the  large-scale  equations  of  ATD.  We  provided  background 
information  and  motivation  for  this  approach,  and  we  discussed  the  basic  approaches  being 
employed  to  correlate  the  chaotic  maps  to  be  used  as  turbulence  models  to  date.  We  then  gave 
some  preliminary  results  for  pipe  flow.  In  the  present  report  we  will  supply  results  from  the 
continuation  of  these  studies.  Considerable  progress  has  been  made,  both  in  terms  of  completed 
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results  and  with  respect  to  analytical  studies  to  put  this  approach  on  a  firm  theoretical  foundation. 


A  final  topic  reported  in  the  Annual  Report  involved  the  ejects  of  Reynolds  averaging. 
These  studies  have  been  concluded  for  the  incompressible  N.-S.  equations,  and  a  papa*  has  been 
submitted  to  the  Physics  of  Fluids  A  (now  Physics  of  Fluids.  ) 

2.  COMPLETE  ADDITIVE  TURBULENT  DECOMPOSITION 

One  of  the  main  problems  proposed  for  study  under  the  current  grant  was  implementation 
of  ATD  for  the  incompressible  N.-S.  equations  in  generalized  coordinates.  The  method  requires 
no  averaging,  or  modeling  at  any  level  due  to  additive  two-scale  decomposition  of  governing 
equations.  Thus,  like  direct  numerical  simulation,  it  is  completely  consistent  with  the  original 
unaveraged  equations;  but  required  arithmetic  is  significantly  reduced  via  consistent  linking  of 
large-scale  and  small-scale  phenomena,  resulting  in  the  ability  to  focus  on  local  regions.  For 
completeness,  we  will  summarize  the  main  idea  of  ATD  in  Sec.  2.1.  Next,  the  progress  on  the 
large-scale  equations,  mainly  focusing  on  filtering  of  the  aliasing  error  which  arises  from  not 
being  able  to  resolve  all  wavenumber  components  in  the  coarse  large-scale  grid,  is  reported  in 
Sec.  2.2.  In  Sec  2.3,  we  summarize  the  extensive  investigation  of  the  small-scale  equations, 
most  of  which  had  already  been  completed  at  the  time  of  the  Annual  Report,  and  in  Sec  2.4  we 
discuss  our  recent  efforts  pertaining  to  construction  of  global  small-scale  solutions  by  applying 
domain  decomposition-like  algorithms  to  the  collection  of  local  smaU-scale  results. 

2. 1  Theoretical  Formulation  of  ATO 

The  equations  we  are  studying  are  the  2-D  incompressible  N.-S.  equations  given  as 


VU  =  0, 


(1) 


4^-»-U-VU=  -VP+— AU. 
dt  Re 


(2) 
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Here  U  =  (U.V)^  is  the  velocity  field,  and  P  is  pressure(  divided  by  constant  density);  A  is  the 
2-D  Laplacian.  Using  additive  decomposition  we  first  decompose  dependent  variables  in  the 
usual  manner: 


U  =  u  +  u*  ,  V  =  V  +  V*  ,  P  =  p  +  p*.  (3) 

where  u,  v  and  p  are  large-scale  quantities,  and  “  *  *’  denotes  the  small-scale  part.  It  should  be 
remarked  that  large-scale  quantities  are  not  averages,  and  should  be  interpreted  as  the  first  few 
terms  of  a  Fourier  representation;  the  small-scale  part  then  corresponds  to  the  series  remainder. 
Substituting  Eqs.  (3)  into  Eq.  (1),  we  get  the  large-  and  small-scale  continuity  equations, 

u»  +  Vy*-»^,  (4a) 


U*  +  V*  s 

»  y 


-5, 


(4b) 


where  S  is  the  decomposition  divergence,  which  here  we  set  to  zero.  We  note  that  the  choice 
5  =  0  is  not  necessary  within  the  ATD  formalism;  but  it  is  convenient,  and  it  allows  ATD  to  more 
nearly  conform  with  classical  approaches.  Applying  ATD,  as  analyzed  by  McDonough  et  al. 
(1989),  we  obtain  the  coupled  system  of  differential  equations  for  u  and  u"*: 


~  +  V*(uu)  +  V-(u*u)  =  -  Vp  -I-  t^Au, 
at  Re 


(5a) 


fill*  1 

+  V-(uu*)  -I-  V-(u*u*)  =  -  Vp*  +  ---Au*. 
ot  Re 


(5b) 


We  note  that  Eqs.  (4,  5)  comprise  six  equations  for  the  six  unknowns;  hence  there  is  no  closure 
problem,  and  no  need  for  modeling.  Equation  (5a)  is  the  equation  for  the  large-scale  velocities, 
and  it  can  be  readily  solved  by  standard  numerical  methods,  as  described  in  the  next  section. 

One  of  the  main  features  of  the  overall  ATD  algorithm  as  originally  proposed  by 
McDonough  et  al.  (1984a,  1984b)  is  the  spatial  decomposition  of  the  small-scale  equation  (5b) 
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into  local  systems  defined  on  non-overlapping  subdomains  containing  each  large-scale 
discretization  point.  This  is  shown  schematically  in  Fig.  1.  The  local  equations  can  be  solved 
independently,  providing  complete  parallelization,  and  recoupled  to  obtain  the  global  small-scale 
solution  via  domain  decomposition  techniques  (cf.  Glowinski  et  al.,  1988).  We  remark  that 
although  there  are  opportunities  for  parallelization  in  most  algorithms  associated  with  the 
standard  methods,  these  usually  lie  at  the  numerical  analytic  level.  However,  ATD  provides  a 
high  degree  of  natural  parallelizability  in  its  basic  structure,  completely  independent  of  speciHc 
numerical  methods  used  in  the  solution  process  (and,  of  course,  the  numerical  methods  sdll  can 
be  parallelized).  Thus,  ATD  will  provide  advantages  on  massive  parallel  processors  (MPPs)  that 
are  not  available  in  standard  approaches  to  turbulence  simulation. 


n 

± 


Sma)l-Seal*  Qalafliln  SuMomaIn 
In  Waignbortiood  of  (IJ)  OfW  Point 


Figure  1.  Local  Galerkin  Region 


2.2  Large-Scale  Calculations 

The  large-scale  equations  have  been  coded  in  the  generalized  coordinates  based  on 
Gresho’s  projection  method  (1990),  The  equations  are  discretized  using  centered  differencing  in 
space  and  generalized  trapezoidal  integration  in  nme.  The  details  of  implementation  can  be 
found  in  the  Annual  Report  by  McDonough  et  al.  (1992). 

It  is  well  known  that  the  reason  for  filtering  the  N.-S.  equations  in  large  eddy  simulation 
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(LES)  is  to  produce  a  system  of  equations  whose  solutions  can  be  completely  resolved  with  a 
reasonable  discretization.  The  price  that  is  paid  for  this  is  introduction  of  an  eddy  viscosity  that 
leads  to  non-physical  dissipation  terms  in  the  differential  equations,  and  the  need  for  subgrid- 
scale  (SGS)  modeling  for  constructing  this  eddy  viscosity.  In  the  ATD  algorithm  we  have 
deliberately  avoided  averaging  or  filtering  the  governing  equations,  but  this  means  that  we  must 
Hnd  a  way  to  deal  with  aliasing  errors  that  arise  from  not  being  able  to  resolve  all  wavenumber 
components  present  in  the  actual  solution. 

As  mentioned  earlier,  the  frnite-difference  method  is  used  to  solve  the  large-scale 
equations,  because  it  is  simple  to  implement  and  also  the  most  commonly  used  method  in  current 
application  codes.  It  is  well  known  that  oscillatory  behavior  may  occur  in  the  centered  finite 
difference  solution  when  the  cell  Reynolds  number  is  greater  than  two  (Roach,  1972).  The  cell 
Reynolds  number  is  defined  as 


where  u  is  local  velocity;  p  is  dynamic  viscosity,  and  h  is  grid  spacing.  This  results  in  the  loss  of 
diagonal  donunance  and  positivity  of  the  associated  tridiagonal  matrices,  and  consequent 
oscillatory  fundamental  solutions  to  the  difference  equations.  The  oscillations  can  be  eliminated 
by  reducing  the  grid  size  and  therefore  the  local  value  of  Re„„.  But  this  is  not  always  possible 
due  to  core  memory  and  CPU  time  requirements.  The  spurious  oscillations  can  also  be 
controlled  by  using  schemes  that  are  intrinsically  dissipative,  such  as  the  upwind  differencing,  or 
by  explicitly  adding  artificial  viscosity  to  damp  out  numerical  oscillation.  However,  upwind 
differencing  alters  the  difference  equations,  and  in  extreme  cases  (such  as  simple  first-order 
upwinding)  results  in  qualitatively  incorrect  solutions;  and  addition  of  artificial  dissipation 
formally  changes  both  the  differential  and  difference  equations. 

The  other  possible  approach  is  to  use  frlters  to  suppress  the  numerical  oscillations,  as  has 
been  successfully  done  in  shock  capturing  (Engquist  et  al.,  1989  and  Shyy  et  al.,  1992).  As 
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suggested  by  Shyy  et  al.  (1992),  one  may  choose  to  extract  the  “useful”,  i.e.,  physically 
realizable,  information  from  oscillatory  solutions  obtained  using  unsatisfactory  numerical 
schemes.  The  idea  is  to  eliminate  undesirable  portions  of  the  solution  while  retaining  only  the 
desired,  physically  realizable  ones.  This  approach  is  fundamentally  different  from  LES  where  a 
filter  is  used  to  filter  the  differential  equations  rather  than  the  «olutions.  As  a  consequence,  there 
is  no  inherent  closure  problem,  and  hence,  no  required  modeling  when  solutions  are  Hltered. 

Motivated  by  the  work  of  Shuman  (1957)  and  Shapiro  (1970),  we  choose  to  use  a  family 
of  Shuman-like  filters  designed  as  a  postprocessor  to  damp  the  information  that  is  not  resolvable 
on  the  large-scale.  The  2-D  filters  are  in  the  form 


+  f . .  ,  +  f . . , 


4  +  m 


+  mf. . 
_ 'A 


(6) 


These  filters  are  second-order  accurate  with  a  truncation  error  (h*u„  +  hyUyy)/(4  +  m).  In  general, 
m  can  take  on  any  value  greater  than  four.  However,  the  optimal  choice  is  governed  by  the 
minimum  amount  of  numerical  viscosity  required  to  suppress  the  aliasing  error  growth.  As 
pointed  out  by  Khosla  and  Rubin  (1980)  there  may  be  no  general  way  of  arriving  at  an  optimal 
value  of  m  analytically;  it  must  be  obtained  by  numerical  experiments.  The  analysis  of  above 
family  of  filters  shows  the  following  properties:  i)  they  provide  good  results  on  reasonably  fine 
grids,  and  ii)  they  have  little  effect  on  the  low  wavelength  components  of  the  solutions. 


We  have  conducted  a  number  of  numerical  experiments  to  aid  in  our  characterization  of 
the  behavior  of  the  filters  represented  by  Eq.  (6).  These  have  included  studies  of  both  steady  and 
time-dependent  1-D  Burgers’  equations  and  in  conjunction  with  time-accurate  projection 
methods  for  the  lid-driven  cavity  problem.  Figure  2  shows  the  filtered  and  unfiltered  solution  for 
a  steady  Burgers’  equation  problem  with  the  exact  solution  given  as 

u  =  ^(1  -  tanh-^). 

2  4v 
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For  these  results,  the  g:id  spacing  and  the  pseudo-time  step  size  are  both  0.2.  The  diffusion 
coefficient  is  taken  to  be  10'^.  Hence,  the  cell-Reynolds  number  restriction  is  grossly  violated. 
The  smooth  solution  is  obtained  using  a  filter  with  m  =  3.4  in  Eq.  (6),  which  is  almost  identical 
to  the  exact  solution.  Without  filtering,  the  solution  diverges.  The  oscillatory  solution  shown  is 
obtained  using  a  filter  with  m  =  100,  which  corresponds  to  a  very  weak  filter.  We  observe  that 
with  an  appropriate  filter,  oscillations  characteristic  of  large  cell  Reynolds  numbers  can  be 
eliminated,  i.e.,  the  aliasing  error  growth  can  be  controlled.  The  fact  that  the  filtered  solution  is 
almost  identical  to  the  exact  solution  demonstrates  that  the  filters  are  able  to  distinguish  the 
aliased  solution  from  the  unaliased  solution. 


Figure  2.  Comparison  of  Solution  Profiles  with  Different  Filters 

The  above  analysis  is  confined  to  the  stationary  solution  of  Burgers’  equation.  The  next 
problem  considered  was  a  time-dependent  Burgers’  equation  with  forcing  term  such  that  the 
exact  solution  becomes  u  =  x‘.  For  short  time,  no  oscillation  is  found  because  local  the  cell 
Reynolds  number  is  small.  When  we  integrate  Burgers’  equation  for  a  longer  time  with  no  filter, 
zigzag-type  oscillation  appears  in  the  solution.  Figure  3  shows  part  of  the  filtered  and  unfiltered 
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Comparison  of  Solution  Profiles  with  and  without  Filtering 


Figure  4.  Filtered  Velocity  Field  of  Lid-Driven  Cavity  at  Re  *  5000 


nunwrical  solution  as  well  as  the  exact  solution  at  1 3:6.0 .  The  solutions  are  obtained  with  both 
grid  spacing  and  time  step  being  O.OOS.  Once  more  it  is  demonstrated  that  filtering  can  suppress 
the  growth  of  errors  arising  from  aliasing. 

Numerical  calculations  were  also  performed  for  the  2-D  lid-driven  square  enclosure  flow 
with  and  without  filters.  Figure  4  presents  velocity  vector  plots  for  Re  =  5000  on  a  41x41  grid. 
Without  filters  the  solution  diverges  due  to  the  aliasing  error  growth.  This  plot  agrees 
qualitatively  with  published  results,  showing  the  primary  central  vortex  and  secondary  vortices 
in  both  of  the  two  lower  comers. 

The  complete  report  on  the  study  of  the  linear  filters  can  be  found  in  a  forthcoming  paper 
by  Yang  and  McDonough  (1994).  We  would  like  to  conclude  the  discussions  presented  here  by 
making  the  following  remarks  concerning  the  numerical  results  given  in  this  section.  First,  quite 
satisfactory  solutions  have  been  obtained  from  otherwise  oscillatory  ones  via  the  filtering 
procedure.  Second,  for  the  steady  state  problems,  we  need  to  add  a  pseudo-time  derivative  in  the 
equation  to  be  solved,  and  use  a  time-stepping  method  to  get  the  steady-state  solution.  This  idea 
is  similar  to  the  multistage  smoothing  in  the  muidgrid  method.  As  implied  by  Majda  et  al. 
(1978),  and  also  pointed  out  by  Shyy  et  al.  (1992),  it  is  too  late  just  to  filter  the  final  steady-stote 
solution  (if  one  is  obtained)  because  the  characteristic  wavelength  of  the  oscillations  is  too  long 
to  be  affected  by  a  low-pass  filter.  From  the  discussion  in  the  previous  section,  we  know  in  this 
case  the  filter  can  do  little  to  improve  the  solution  accuracy.  On  the  other  hand,  the  high- 
frequency,  short  wavelength  oscillations  may  lead  to  nonlinear  instabilities,  and  divergence, 
before  a  steady  state  can  be  achieved  if  the  filter  is  not  applied  at  every  time  step.  Finally,  we 
should  note  that  the  use  of  post-processing  filters  is  particularly  attractive  in  the  context  of  the 
ATD  algorithm  because  if  the  filter  can  be  accurately  characterized  so  that  it  is  known  how  much 
of  the  large-scale  information  has  been  removed,  it  is  possible  to  add  this  back  in  during  the 
small-scale  calculations.  Thus,  it  may  be  of  value  to  employ  ATD-like  algorithms  even  for 
laminar  flow  calculations  as  a  means  of  treating  cell  Reynolds  number  problems. 
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We  have  made  steady  progress  toward  completion  of  the  large-scale  calculation  procedure. 
The  current  code  is  undergoing  thorough  investigation  for  lid-driven  cavity  flow  in  generalized 
coordinates,  and  generation  of  the  grid  for  flow  over  a  NACA  0012  airfoil  is  finished.  The  proper 
implementation  of  outflow,  inflow  and  solid  boundary  conditions  is  now  in  progress,  and  more 
efficient  Poisson  solvers  based  on  the  strongly  implicit  procedure  (SIP)  are  being  implemented 
and  compared  with  SOR  and  ADI.  We  expect  this  will  greatly  reduce  the  CPU  time  required  for 
each  time  step  since  most  CPU  time  is  spent  on  solving  the  Poisson  equation  resulting  from  use 
of  projection  methods. 

2.3  Small-Scale  Calculations 

The  complete  investigation  of  small-scale  performance  and  the  large-  to  small-scale 
transfers  of  data  is  reported  in  Yang  and  McDonough  (1992a).  Unlike  the  large-scale  calculation, 
the  small-scale  equations  are  solved  through  constructing  local  Galerkin  approximations  in  the 
neighborhood  of  each  finite  difference  grid  point,  as  indicated  in  Fig.  1. 

We  assume  the  solution  can  be  represented  by  the  truncated  Fourier  series  defined  as 
follows  on  a  local  domain  containing  the  grid  point  (Xi,yp; 

K 

u*(x,y,t)  =  ^  j^\t)cosa^x*sina^y*,  (7a) 

lun=l 

K 

vj(x,y,t)=  ^b[^\t)sina^x*cosa^y*,  (7b) 

kjn=l 

K 

pj(x,y,t)=  ^c|;^\t)sina^x*sina^y*,  (7c) 

kjn=I 


where 
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X*  =  X  -  (x.  -  h/2),  y*  =  y-(y--  h/2). 

with  (x.y)  €  (Xi-h/2,  Xj+h/21x(y.-h/2,  y.+h/2],  and  a„  =  nui/h.  Here  h  is  the  finite  difference 
grid  spacing  for  the  large-scale  calculation. 

By  substituting  Eqs.  (7)  into  (4b)  and  (Sb)  and  constructing  the  Galerkin  inner  products,  we 
obtain  a  system  of  ordinary  differential  equations  for  the  time-dependent  Fourier  coefficients. 
Details  can  be  found  in  Yang  and  McDonough  (1992a,b).  The  local  large-scale  velocity  data  are 
directly  input  to  the  small-scale  equations  and  appear  as  coefficients  of  the  convolution  terms. 
This  allows  computation  of  the  small-scale  solution  without  explicitly  representing  all  modes  of 
the  large-scale  as  well.  This  is  the  critical  difference  between  small-scale  ATD  and  usual  DNS, 
and  it  renders  ATD  far  more  efficient  for  local  calculations  than  would  be  DNS.  At  the  same 
time  we  believe  it  is  capable  of  providing  more  faithful  simulations  than  does  “minimal”  DNS 
proposed  by  Moin  (e.g.,  Moin,1992),  specifically  because  of  this  large-scale  information. 

The  resulting  system  of  ordinary  differential  equations  are  integrated  using  a  second-order 
Runge-Kutta  method  (Heun’s  method)  with  constant  time  step.  Complete  details  of  choice  of 
initial  conditions,  etc.  can  be  found  in  Yang  and  McDonough  (1992a,b).  Figure  5  displays  a 
typical  time  series  corresponding  to  a  chaotic  solution.  The  value  of  Re  at  which  this  calculation 
has  been  performed  (Re  =  10^)  is  several  orders  of  magnitude  higher  than  has  previously  been 
possible  with  standard  methods.  This  is  accomplished  by  performing  the  calculations  on  very 
small  length  scales,  and  by  supplying  large-scale  input  to  these  scales  in  a  consistent  manner. 
Nevertheless,  it  is  important  to  check  the  resolution  adequacy  of  these  results.  Figure  6  provides 
an  indication  of  the  degree  of  convergence  of  the  small-scale  solution  by  depicting  the  rate  of 
decay  of  Fourier  coefficient  time-averaged  amplitudes  as  a  function  of  wavenumber.  In  addition 
we  comment  here  that  a  comparison  of  our  overall  effective  resolution  (large-scale,  based  on  the 
small-scale  subinterval  length,  plus  small-scale  modal  representation)  with  predictions  by  Foias 
and  Treve  (1981)  regarding  required  number  of  modes  as  a  function  of  Reynolds  number  shows 
that  our  calculations  should  be  fairly  reliable. 
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Figure  6. 


Our  experience  working  with  the  small-scale  equations  has  shown  that  the  discrete 
convolution  evaluations  consume  a  tremendous  fraction  of  computing  time,  as  would  be 
expected,  since  this  requires  0(K*)  arithmetic  per  time  step  for  each  of  the  K  ODEs.  We  also 
note  that  the  number  of  modes  typically  used  in  the  small-scale  calculation  is  not  large  enough 
for  a  FFT  to  provide  much  improvement  The  remedy  is  to  take  advantage  of  rapid  development 
of  parallel  computers.  Our  preliminary  parallel  processing  studies  using  6-processor  IBM 
3090-600JS  demonstrate  that  speedup  is  linear  with  number  of  available  processors.  At  the  same 
time,  we  expect  to  achieve  significant  but  less  than  linear,  speedup  by  parallelizing  ODE 
solvers  on  each  small-scale  local  domain.  This  woric  is  scheduled  to  start  very  soon.  Fulfrllment 
of  this  task  will  enable  us  to  simulate  turbulent  flows  at  high  Reynolds  numbers  that  are  not 
achievable  by  LES  and  DNS. 

2.4  Domain  Decomposition 

We  have  as  yet  not  implemented  a  specific  domain  decomposition  procedure  for  the  small- 
scale  equations,  but  we  have  been  studying  two  different  approaches.  Results  due  to  Chan  and 
Mathew,  reported  in  McDonough  et  al.  (1989)  suggest  that  it  may  not  be  necessary  to  employ  a 
traditional  domain  decomposition  algorithm  in  the  context  of  ATD  because  the  small-scale 
solutions  are  needed  only  at  large-scale  grid  points,  which  happen  to  lie  well  within  the  interior 
of  the  small-scale  subdomains,  and  because  the  effects  of  boundary  errors  (the  thing  that  must  be 
corrected  by  domain  decomposition  iterations)  are  not  very  large  in  the  interior  of  the 
subdomains. 

Our  first  approach  to  recoupling  the  local  small-scale  calculations  to  obtain  the  global 
small-scale  solution  makes  direct  use  of  this.  This  method  is,  in  a  sense,  a  predictor-corrector 
procedure  requiring  two  separate  solutions  of  the  small-scale  equations.  The  prediction  is  made 
using  overlapping  subdomains  of  size  2h,  where  h  is  the  large-scale  grid  spacing,  and  Fourier 
representations  containing  K+1  modes.  These  calculations  can  be  performed  in  an 
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“embarrassingly”  parallel  fashion  with  no  iteration,  and  there  are  essentially  the  same  number  of 
such  subdomains  (in  1-D)  as  occur  for  the  fully-resolved  h-subdomain  small-scale  calculations. 
It  can  be  seen  from  Fig.  7  that  the  center  of  each  of  the  2h-subdomains  (where  errors  are 
smallest)  coincides  with  the  boundary  between  two  h-subdomains.  Thus,  the  2h-subdomain 
calculation  should  provide  a  reasonably  accurate  prediction  of  boundary  values  for  the  h- 
subdomain  calculations.  The  latter  are  then  performed,  again  non-iteratively,  and  their 
results — taken  from  the  center  of  the  h-subdomains — are  used  to  construct  the  small-scale 
solutions  for  use  in  the  large-scale  equations.  It  should  be  observed  that  at  each  stage  of  this 
process,  results  that  are  actually  used  are  taken  from  the  part  of  the  subinterval  at  which  errors 
should  be  a  minimum.  This  predictor-corrector  form  of  domain  decomposition  is  currently  being 
tested  in  the  context  of  1-D  model  problems. 

2h-Grid  Overlapping  Smatl-Scala  Subdomaina 


U 


Figure  7.  Domain  Decomposition  for  Recoupling  Small-Scale  Solutions 

The  second  method  that  we  are  considering  for  recoupling  the  small-scale  solutions  makes 
use  of  the  well-known  fact  from  domain  decomposition  theory  that  the  fine-grid  recoupling 
iterations  converge  faster  if  there  is  an  underlying  coarse  grid  calculation.  (This  fact  is  used 
implicitly  in  the  first  approach  also.)  In  particular,  after  the  small-scale  solutions  have  been 
computed  on  all  of  the  individual  subdomains  independently,  they  can  be  combined  to  form  a 
global  solution  in  the  following  way.  We  discretize  the  small-scale  equations  (5b)  on  the  large- 
scale  time  step.  Let  n-t-1  be  the  desired  advanced  large-scale  time  level,  and  suppose  N*  is  the 
small-scale  time  step  corresponding  to  this.  Let 
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„-'  =  „"+i[^(„-0  +  L,(u’)] 

be  a  generic  representation  of  Eq.  (Sb).  Then  in  terms  of  the  above  notation  we  have 

where  k,  is  the  small-scale  time  step,  and  h  is  the  large-scale  space  step.  We  observe  that  Eq.  (8) 
directly  leads  to  nearest-neighbor  coupling  of  the  small-scale  solutions,  but  in  the  context  of 

*  N* 

projection  methods  should  probably  be  preceded  by  a  mass-conservation  correction  of  both  u  ' 

and  u '  ' .  Thus,  completely  global  coupling  will  be  achieved  through  the  pressure  Poisson 
equations.  In  particular,  Eq.  (8)  could  be  used  iteratively  to  set  boundary  conditions  for  the  next 
set  of  small-scale  calculations,  but  it  seems  unlikely  that  this  would  be  necessary. 

It  is  this  second  procedure  that  we  plan  to  employ  first  in  the  construction  of  our  complete 
ATD  algorithm  due  to  its  relative  simplicity  for  muld-dimensional  calculations.  We  also  note  in 
closing  this  section  that  for  high-Re,  fully  developed  turbulence,  it  may  be  unnecessary  to 
expend  much  effort  in  recoupling  the  small-scale  solutions  because  spatial  correlation  lengths 
may  be  far  smaller  than  the  large-scale  grid  spacings. 

The  final  theoretical  area  that  has  received  attention  by  the  PI  concerns  using  very  local 
(only  a  few  grid  points  in  each  direction)  finite  difference  (or  finite  volume)  discretization  of  the 
small-scale  equations.  These  are,  of  course,  nonlinear  algebraic  maps,  but  not  a  great  deal  is  yet 
known  regarding  their  properties.  Such  a  formalism  lies  between  the  large-scale  ATD/chaotic 
map  model  approach  and  complete  ATD.  It  is  clearly  consistent  with  the  N.-S.  equations,  and  it 
is  embarrassingly  parallel.  At  the  same  time,  it  is  far  more  efficient  than  complete  ATD  (but 
would  converge  to  complete  ATD  as  the  small-scale  domain  sizes  approach  the  large-scale  grid 
spacing,  or  vice  versa),  and  unlike  complete  ATD  it  would  be  fairly  straightforward  to 
implement  in  the  context  of  adaptive,  multi-level  (more  than  two  levels)  unstructured  grids.  We 
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have  so  far  been  able  to  demonstrate  only  what  would  be  expected  on  intuitive  grounds;  namely, 
the  accuracy  of  the  large-scale  grid  spacing  can  be  maintained,  and  qualitatively  realistic  small- 
scale  solutions  can  be  computed  to  this  same  level  of  accuracy. 

3.  LARGE-SCALE  ATD  WITH  CHAOTIC  MAP  TURBULENCE  MODELS 

In  this  section  we  will  discuss  an  aspect  of  ATD  that  was  new  to  the  current  grant,  namely 
use  of  the  large-scale  equations  (4a,  Sa)  alone  with  fluctuating  quantities  modeled  with  chaotic 
algebraic  maps.  We  have  provided  considerable  information  on  this  approach  in  the  1992  Annual 
Report  (McDonough  et  al.,  1992),  so  the  present  discussion  will  focus  mainly  on  results  obtained 
during  1993.  Nevertheless,  we  will  include  some  basic  background  information  in  this  report  for 
the  sake  of  completeness.  Following  this  background  information  will  be  a  subsection  in  which 
we  present  further  results  obtained  by  applying  this  approach  to  pipe  flow.  Then  in  a  final 
subsection  we  will  discuss  some  theoretical  developments  associated  with  method. 

3.1  Background 

We  first  remark  that  in  light  of  recent  results  showing  that  Reynolds-averaged  approaches 
cannot  be  consistent  with  the  Navier-Stokes  equations  (McDonough,  1993),  and  that  this  lack  of 
consistency  results  from  having  to  model  terms  that  arose,  in  the  first  place,  as  a  consequence  of 
averaging,  it  would  seem  reasonable  to  investigate  the  possibility  of  developing  turbulence 
models  based  on  unaveraged  equations.  ATD  provides  a  quite  natural  setting  in  which  this  can  be 
done.  In  particular,  we  can  see  from  Eq.  (Sa)  that  if  u*  can  be  modeled  in  such  a  way  as  to 
approach  zero  as  the  grid  spacing  goes  to  zero,  then  this  approach  would  be  consistent  with  the 
N.-S.  equations.  It  is  interesting  to  note  that  the  Reynolds  stresses  in  LES  typically  have  this 
property,  but  there  remains  a  realizability  problem  in  LES,  at  least  when  dynamic  SGS  models 
are  used,  that  is  completely  absent  in  ATD.  In  fact.  Galilean  invariance  and  realizability  are  all 
but  automatic  in  ATD  because  it  is  the  fluctuating  quantities  themselves  that  are  modeled.  This  is 
a  major  advantage  of  this  approach. 
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3.2  Intennittent  Turbulent  Pipe  Flow 


As  we  have  reported  in  the  1992  Annual  Repoit,  we  have  been  studying  large-scale  ATD 
with  chaotic  map  models  for  intennittent  turbulence  in  an  axisymmetiic  pipe.  Although  it  is  not 
entirely  realistic  for  turbulent  flow,  we  have  retained  the  axial  symmetry  in  our  formulation.  The 
basic  algorithm  proceeds  in  the  following  way  for  each  time  step: 

1.  Compute  an  estimate  of  u  from  (4a,  Sa)  using  u*  brom  time  level  n. 

2.  Use  these  large-scale  results  as  parameters  in  chaotic  maps  to  calculate  at  each  grid 
point. 

3.  Correct  u"*'  by  solving  (4a,  5a)  with  in  the  crossterms. 

4  Construct  the  complete  solution  at  time  level  n-f  1: 


n+l  n+l  .  B  +  1  n  +  1  ,  *J>+1 

U  =!U  +u  ,  p  =p  +p  . 

It  should  be  observed  that  in  our  current  algorithm  we  apply  a  mass  conservation  calculation  to 
the  small-scale  velocity  fluctuations  created  from  the  algebraic  map  in  order  to  guarantee  that 
V  u*  =  0  (and  thus  V  U  =  0),  and  to  simultaneously  generate  p’’"*‘  for  calculating  p"*‘.  We  do  not, 
however,  guarantee  that  the  small-scale  momentum  equations  (Sb)  are  satisfied,  either  locally  or 
globally.  We  have  discussed  in  McDonough  et  al.  (1992)  the  various  criteria  we  feel  must  be  met 
to  assure  that  the  chaotic  map  models  behave  in  a  realistic  way,  but  even  if  all  such  criteria  are 
satisfied  we  cannot  be  guaranteed  that  u*  is  truly  close  to  a  solution  to  (Sb).  Thus,  the  procedure 
as  we  have  implemented  it  to  date  is  indeed  a  model,  and  although  it  is  guaranteed  to  be 
consistent  with  the  N.-S.  equations  as  discretization  step  sizes  approach  zero,  little  can  be 
guaranteed  a  priori  for  the  large  step  sizes  one  might  typically  use  in  solving  practical 
engineering  problems. 


From  a  physical  viewpoint  we  see  that  Eq.  (5a)  accounts  for  transport  of  small-scale 
quantities  by  the  large-scale  velocity  field,  but  not  vice  versa.  It  is  possible  that  in  many  flow 
situations  this  is  reasonably  accurate,  and  for  such  flows  this  modeling  procedure  may  produce 
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quite  accurate  results.  But  in  the  contrary  case  it  ntay  not  It  should  already  be  clear,  however, 
that  there  is  an  easily  implemented,  relatively  inexpensive  improvement  that  can  be  made  to 
cover  this  problem.  Namely,  construct  a  global  small-scale  solution  in  the  same  manner  as 
discussed  in  Sec.  2.4,  but  now  start  with  results  from  chaotic  maps  rather  than  local  solutions  to 
the  N.-S.  equations.  We  have,  as  yet,  not  implemented  this  promising  addition,  so  the  results  we 
now  present  were  obtained  from  our  original  formulation. 


Figure  8.  Velocity  Fluctuation  Distribution  from  Chaotic  Map  Model 
for  Iteration  step  =  1500, 3000 

These  calculations  were  performed  on  a  quite  coarse  41x21  grid  for  a  pipe  with  L/D  =  20 
and  Re  =  2000.  Thus,  it  can  be  expected  that  the  flow  will  not  be  fully  developed  at  the 
downstream  end  of  the  pipe.  With  the  uniform  grid  spacing  used  in  the  radial  direction,  we 
expect  to  have  no  more  than  one  grid  point  in  the  log  layer  near  the  pipe  wall,  but  the  goal  in  the 
present  study  is  more  to  explore  the  potential  of  the  method  rather  than  to  produce  highly 
accurate  results.  Figure  8  displays  the  spatial  distribution  of  fluctuating  velocities  at  two  different 
times  during  the  calculations.  The  chaotic  map  being  employed  here  is  simple  the  logistics  map. 
It  should  be  noted  that  the  velocity  vectors  are  plotted  on  a  much  larger  scale  than  that  used  in 
the  calculations.  In  particular,  u*  -  0.  lu  is  rather  typical.  Close  examination  of  Fig.  8  reveals  the 
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spatial  and  temporal  chaos  that  would  be  expected  in  a  turtnilent  flow,  as  well  as  indications  of 
intermittency. 

Figure  9  displays  the  distribution  of  large-scale  velocities  at  the  same  times  as  in  Fig.  8. 
We  note  that  if  the  fluctuating  velocities  were  identically  zero  ,  the  large-scale  velocities  should 
be  very  regular  (although  not  fully  developed  for  Re  =  2(XX)).  But  die  effects  of  the  chaotic 
fluctuations  are  clearly  evident.  Finally,  in  Fig.  10  we  present  time  averaged  velocity  profile 
along  the  pipe.  This  has  been  obtained  by  time  averaging  the  results  computed  for  the  complete 
velocities  at  every  time  step.  It  can  be  seen  that  the  flow  is  still  developing,  and  that  the  velocity 
profile  is  distinctly  different  from  the  laminar  case.  However,  the  distribution  of  turbulent  kinetic 
energy  is  not  in  complete  agreement  with  what  is  expected  for  turbulent  pipe  flow,  so  evidently 
additional  work  needs  to  be  put  into  development  of  the  chaotic  maps.  More  details  regarding 
these  calculations  will  appear  in  a  forthcoming  paper,  McDonough  and  Zhong  (1994). 

3.3  Theoretical  Studies 

It  is  not  difficult  to  see  that  there  are  numerous  theoretical  aspects  of  the  large-scale 
ATD/chaotic  map  modeling  approach  that  should  be  dealt  with.  We  will  mention  a  few  here. 

First,  we  would  expect  on  intuitive  grounds  that  the  amplitudes  of  the  small-scale 
fluctuations  should  depend,  in  some  way,  on  local  large-scale  flow  properties,  but  the  question  of 
how  to  quantify  this  dependence  needs  to  be  addressed.  Moreover,  as  we  have  already  indicated, 
we  must  require  that  these  amplitudes  also  depend  on  the  (local)  discretization  step  sizes,  and  in 
particular  approach  zero  with  the  discretization  step  size.  But  these  are  all  rather  vague, 
qualitative  requirements.  Work  has  recently  been  underway  by  the  PI  and  students  funded  by 
AFOSR  Grant  #49620-92-1-0488  to  provide  quantitative  results  along  these  lines.  Details  will  be 
reported  in  Hylin  et  al.  (1994)  and  Weatherly  et  al.  (1994).  These  studies  are  quite  general  and 
are  being  done  in  the  context  of  2-  and  3-D  compressible  N.-S.  equations.  Some  of  the  results, 
however,  are  quite  generic  and  apply  also  to  incompressible  flows.  In  particular,  we  have  shown 
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Figure  9.  Large-Scale  Velocity  Distribution  for  Turbulent  Pipe  Flow 
for  Iteration  step  =  1500, 3000 


Figure  10.  Averaged  Total  Velocity  Distribution  for  Turbulent  Pipe  Flow 

that  velocity  amplitudes  obtained  from  chaodc  maps  should  scale  according  to  the  following 
formula; 

\  =  h  llVull 

where  h  is  the  large-scale  grid  spacing;  v  is  kinematic  viscosity,  and  is  a  constant  that  must  be 
determined  from  data  for  particular  classes  of  flows. 

A  second  area  of  study  involves  details  of  the  structure  of  the  algorithm  for  implementing 
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large-scale  ATD  with  chaotic  maps.  There  are  various  alternatives  to  the  one  provkled  in  the 
preceding  subsection,  distinguished  primarily  by  the  details  of  maintaining  the  divergence-free 
constraint.  As  noted  earlier,  it  is  not  necessary  that  Uiis  be  satisfied  on  either  the  large-  or  small- 
scale  separately  (but  if  it  is  satisfied  on  one,  it  must  be  on  both),  and  this  suggests  a  potentially 
more  efficient  procedure  since  most  of  the  computational  work  on  both  scales  is  devoted  to 
maintaining  the  divergence-free  condition.  As  an  alternative,  we  might  separately  calculate  u  and 
u*  without  mass  conservation,  and  then  find  P  such  that  V-U  =  0.  Once  P  is  known,  the  complete 
equations  can  be  solved  for  U.  This  would  eliminate  half  of  the  work  involving  pressure  Poisson 
equation  solves. 

4.  SUMMARY 

We  begin  this  summary  section  by  observing  that  although  we  have  not  quite  attained  all 
of  the  goals  set  forth  in  the  proposal  for  this  grant,  we  have  nevertheless  advanced  the 
development  of  ATD  quite  significantly  in  both  of  the  two  main  areas  under  investigation:  i) 
complete  (large-scale  plus  small-scale  ATD)  in  2-D  generalized  coordinates  for  the  N.-S. 
equations,  and  ii)  large-scale  ATD/chaodc  map  turbulence  models.  The  list  of  major 
accomplishments  include  the  following. 

1.  Completion  of  Cartesian  coordinate  large-scale  equation  solution  algorithm,  and  near 
completion  in  generalized  coordinates. 

2.  Thorough  analysis  of  post-processing  filtering  techniques  for  treating  the  cell-Re  problem 
arising  in  coarse-grid  large-scale  calculations,  and  demonstration  that  this  is  a  very  appropriate 
approach  in  the  context  of  ATD 

3.  Completion  of  the  generalized  coordinate  local  Galerkin  small-scale  algorithm,  including 
initial  studies  of  parallelization. 

4.  Analysis  of  the  domain  decomposition  techniques  required  to  obtain  global  small-scale 
solutions. 

5.  Setting  criteria  and  procedures  to  be  used  in  developing  chaotic  map  small-scale  turbulence 
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models. 


6.  Construction  and  testing  of  code  for  modeling  intermittent  axisymmetric  pipe  flow. 

7.  Theoretical  analysis  of  the  small-scale  turbulence  models. 


The  results  of  these  accomplishments  have  been,  or  will  be,  reported  in  the  list  of 
publications,  conference  presentations  and  invited  talks  given  below. 

Conference  and  Journal  Papers 

Hylin,  E.  C.,  McDonough,  J.  M.  and  Weatherly,  D.  C.,  1994,  “Modeling  the  Subgrid  Scale  Flow 
with  a  Chaotic  Map,”  to  be  submitted  to  J.  Fluid  Mech. 

McDonough,  J.  M.,  1993,  “On  the  Effects  of  Modeling  Errors  in  Turbulence  Qosures  for 
Reynolds-Averaged  Navier-Stokes  Equations,”  submitted  to  Phys.  Fluids  A. 

McDonough,  J.  M.  and  Zhong,  X.,  1994,  “Intermittent  Pipe  Flow  Computation  Via  Additive 
Decomposition  of  the  Navier-Stokes  Equations  with  Chaotic  Map  Models,”  submitted  to  for 
ASME  Winter  Annual  Meeting,  1994. 

Weatherly,  D.  C.,  Hylin,  E.  C.  and  McDonough,  J.  M.,  1994,  “Additive  Turbulent 
Decomposition  with  Subgrid  Scale  Chaotic  Maps  for  Compressible  Turbulence  Simulation,”  to 
be  submitted  to  J.  Fluid  Mech. 

Yang,  Y.  and  McDonough,  J.  M.,  1992a,  “Bifurcation  Studies  of  Navier-Stokes  Equations  via 
Additive  Turbulent  Decomposition,”  in  Bifurcation  Phenomena  and  Chaos  in  Thermal 
Convection,  Bau  et  al.  (eds.),  HTD-Vol.  214,  ASME,  New  York. 

Yang,  Y.  and  McDonough,  J.  M.,  1994,  “Studies  of  Linear  Filters  For  Treating  Cell-Re  Problems 
in  Finite  Difference  Schemes,”  to  be  submitted  to  Int.  J.  Comput.  Fluid  Dynamics. 

Conference  Presentations  and  Invited  Talk 

Hylin,  E.  C.,  McDonough,  J.  M.  and  Weatherly,  D.  C.,  1993,  “Modeling  the  Subgrid  Scale  Flow 
with  a  Chaotic  Map,”  Bull.  Amer.  Phys.  Soc.  38, 2304 

McDonough,  J.  M.,  1991,  “Unaveraged  Turbulence  Models  Based  on  the  Large-Scale  Equations 
of  Additive  Turbulent  Decomposition,”  Bull.  Amer.  Phys.  Soc.  36, 2648 
McDonough,  J.  M.,  1992,  “Unaveraged  Turbulence  Models  Based  on  the  Large-Scale  Equations 
of  Additive  Turbulent  Decomposition,”  presented  at  Univ.  of  Southern  Calif.  Aerospace 
Enginering  Seminar,  Los  Angeles,  CA,  Mar.  23. 

McDonough,  J.  M.,  1993,  “Intrinsic  Errors  in  Integrations  of  the  Reynolds-Averaged  Navier- 
Stokes  Equations,”  Bull.  Amer.  Phys.  Soc.  38, 2304 

McDonough,  J.  M.,  Zhong,  X.  and  Xiang,  L.,  1992,  ‘Turbulence  Models  Constructed  from 
Unaveraged  Equations  with  Chaotic  Map  Closures,”  presented  at  SIAM  40th  Annual  Mtg.,  Los 
Angeles,  CA,  July  19-24. 

Weatherly,  D.  C.,  Hylin,  E.  C.  and  McDonough,  J.  M.,  1993,  “Additive  Turbulent 
Decomposition  with  Subgrid  Scale  Chaotic  Maps  for  Compressible  Turbulence  Simulation,” 
Bull.  Amer.  Phys.  Soc.  38,  2266 
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